Surface growth effects on reactive capillary-driven flow: Lattice Boltzmann investigation 
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The Washburn law has always played a critical role for ceramics. In the microscale, surface forces take 
over volume forces and the phenomenon of spontaneous infiltration in narrow interstices becomes of particular 
relevance. The Lattice Boltzmann method is applied in order to ascertain the role of surface reaction and 
subsequent deformation of a single capillary in 2D for the linear Washburn behavior. The aim of the proposed 
investigation is to portrait the complex problem of reactive infiltration of molten silicon into carbon preforms. 
Several geometric characteristics for the capillaries are considered, as well as different infiltration and reaction 
conditions. The main result of our work is that the phenomenon of pore closure can be regarded as independent 
of the infiltration velocity, and in turn a number of other parameters. The instrumental conclusion drawn from 
our simulations is that short pores with wide openings and a round-shaped morphology near the throats represent 
the optimal configuration for the underlying structure of the porous preform in order to achieve faster infiltration. 
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I. INTRODUCTION 



II. COMPUTATIONAL MODELS 



The infiltration of porous carbon (C) by molten silicon (Si) 
is a widespread industrial technique for the fabrication of 
composite materials with enhanced thermal and mechanical 
resistance at high temperatures [1 1. Central to this process are 
of course the wetting properties of the substrate against molten 
Si. Importantly, it appears that capillary forces are strongly af- 
fected by the reactivity at high temperatures between Si and 
C to form silicon carbide (SiC). A large body of research has 
been devoted to the subject, definitely defying conventional 
modeling l2]-21 and experimental Q approaches. Here we 
deal with the ability of the Lattice Boltzmann (LB) method 
[f6 , 7 1 to address the dynamic behavior of capillary infiltration 
coupled to surface growth. Of special interest is the retarda- 
tion induced to the Washburn law by the dynamic reduction 
in pore radius. More generally, the present work can help to 
understand experimental results carried out with Si alloys 1 8 - 
[TU . Precisely, the focus is on the linear Washburn law for a 
single capillary in 2D lfT2l [T3l . This behavior is typical for 
the infiltration of pure Si in carbon preforms |14|. Further- 
more, the Washburn analysis for bulk structures is based on 
the results for a single capillary ifTSl IT6l . In the LB model, 
the surface reaction is treated as a precipitation process from 
supersaturation |9, 11, 17-22|. Our simulations indicate that 
the thickening of the surface behind the contact line actually 
limits the infiltration process. This phenomenon is especially 
marked in the vicinity of the onset of pore closure. For this 
reason, we argue that the primary aspect retarding capillary 
infiltration occurs when surface growth inside narrow inter- 
stices obstructs the flow (here the wetting transition at the con- 
tact line is neglected flO' TT, '231). Importantly, it arises that 
higher infiltration velocities, achieved for example by tuning 
some parameters defining the capillary geometry, have little 
to no influence on the process of pore closure by surface re- 
action. It follows that the distribution of pore size is a critical 
parameter for reactive capillary infiltration, besides other as- 
pects to be discussed in the sequel. 



The distinctive feature of the LB method resides is the dis- 
cretization of the velocity space: only a finite number of ele- 
mentary velocity modes is available. Locally, the movement 
of the fluid particles along a particular direction is accounted 
for statistically by means of distribution functions, in terms of 
which the intervening physical quantities are expressed. Hy- 
drodynamic behavior is recovered in the limit of small Mach 
numbers, coiTesponding to the incompressible limit ll24l . It 
has been shown that the problem of capillary flow is better de- 
scribed by a multicomponent system 1.25-27] , since allowing 
to reduce the evaporation-condensation effect lfT2l[T3ll28U30l . 
On a d2q9 lattice, the evolution of every single fluid compo- 
nent is obtained by iterating the BGK equation lISTI 



f^{r + eiAt,t+At) = 

frir,t) f [fnr,t)-fr\pa,ul^)] / = O- S 



(1) 



The lattice spacing is denoted by Ax and the time increment 
by At; hereafter, without loosing generality Av = = 1 (in 
model units). The ff's are the distribution functions for the 
velocity modes e,. The superscript a — 1,2 designates the 
substance. The discrete velocity modes are defined as follows 

r (0,0) /=o 

ei=l (cos[(/-l)7i;/2],sin[(/-l)7r/2]) i=l-4 
[ \/2(cos[(2/-9)7i:/4],sin[(2/-9)7i:/4]) / = 5-8. 

The equilibrium distribution functions read 



1 +3e/ • u^^ 



! = 0-8 



The coefficients w; are defined as wq = 4/9, w, = 1/9 for / = 
1 - 4 and Wi = 1 /36 for / = 5 - 8 1^ [33] . The local density. 
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Per, and the equilibrium velocity, u^, have the explicit form 



Pa{r.t) = l^fr{r,t) 

1=0 



Tg Fa{r,t) 
Pa{r,t) 



F(, is the resultant force experienced by the fluid component 
cr. The velocity m' is given by 

where we have introduced the d-th component fluid velocity 

1 



Pa{r,t) 



t e: 



The velocity u' satisfies the requirement of momentum con- 
servation in the absence of external forces EsI - ETI . In Eq. [T] 
Tcr is the relaxation time, determining the kinematic viscosity 
V(j = (2T(j — 1)/6. The dynamic viscosity is defined as ||T2| 

l^ = pv = Y,Ha = Y.P<y^<y ■ 

a a 

The total density of the whole fluid is of course p — La Per- 
Collisions between the liquid phase and the solid boundaries 
are treated according to the common bounce-back rule, im- 
plementing the no-slip condition for rough surfaces. Cohesive 
forces in the liquid phase are evaluated by means of the for- 
mula OH 

8 

Fc,air,t) = -GcPa{r,t)Y^WiPa{r + eiAt,t)ei ; 

1=1 

a being the other fluid component with respect to a. The 
parameter Gq determines the interaction strength. Adhesive 
forces between the liquid and solid phases are computed using 
the formula (Ml 



Fads,a{r,t) = -Gads.aPa{r,t)Y,Wis{r + eiAt,t)ei . 

1=1 

The interaction strength is adjusted via the parameter Gads. a- 
The function s takes on the value 1 if the velocity mode e, 
points to a solid node, otherwise it vanishes. 

The LB scheme can encompass convection-diffusion sys- 
tems readily under the assumption of low concentrations. This 
condition implies that fluid flow is not affected by molecular 
diffusion. Solute transport is described by introducing the dis- 
tribution functions gi. For solute transport we use the d2q4 
lattice (2TJ, generally applied also for thermal fields I35l[36l . 
The reduced number of degrees of freedom does not restrain 
the predictive capacity of the model ll37l . As for fluid flow, 
the evolution follows the BGK equation ifTSl - ESl 

gi{r + eiAt,t+At)=gi{r,t)-^[gi{r,t)~g.'^{C,u)] 

i= 1-4. 



The relaxation time for solute transport is denoted Ts, defin- 
ing the diffusion coefficient D — (2Ts — l)/4 [21 J. For the 
d2q4 lattice, the equilibrium distribution functions take the 
form 121] 

g^''(C,M)-^C[l+2e,.w] i^l-4. (2) 

The solute concentration is given by C = Ligi- The local sol- 
vent velocity u is assumed to be the overall fluid velocity ll25l - 



u{r,t)^(Y, Pa{r,t)ug{r,t) 



It may be useful to have interfaces also for solute transport 
in order to prevent diffusion from the liquid to the vapor phase. 
In analogy with multiphase fluid flow 12511381 . we add to the 
velocity u in Eq. |2]a term of the form TsFs/C. The function 

is defined as 

4 

Fs{r,t) = -GMr,t)Y.(p{r + eiAt,t)ei . (3) 

i=l 

The parameter Gs allows to tune the intensity. The function 
(p{r,t) is given by 

(p(r,f) = (Poe"^°/^^''''' . 

This mechanism introduces an interface acting as a barrier to 
diffusion moving with the overall fluid. Far from interfaces. 
Fa vanishes and diffusive transport is restored. The contribu- 
tion of i^s is taken into account also for the nodes belonging 
to the surface of the solid phase. 

At the reactive boundary, the solute concentration C satis- 
fies the macroscopic equation llT8H22ll 

an ' 

where is the reaction-rate constant and Cs the saturated 
concentration. From the constraint of the above equation, it 
is possible to calculate the solute concentration at the solid- 
liquid interface 11211 



C 



2gout+feiGs 



where gout is the distribution function of the modes leav- 
ing the fluid phase. At the solid boundaries, the distribution 
functions g\n for the modes e, pointing in the fluid phase re- 
main undetermined. They are updated according to the rule 

gin — gout 

+ C/2; gin and gout are associated with opposite 
directions. As time goes on, solute mass deposits on the react- 
ing solid surface. The mass is updated iterating the equation 

b{r,t+At)^b{r,t)+k,{C~C^) . 
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Figure 1: Scaled density pa/po as a function of GcPo- Tlie upper 
branch corresponds to Pi/po inside the bubble (main density), while 
the lower branch is pa/po (dissolved density). Inset: Dependence of 
the scaled surface tension y/po on GqPq. The results of the simula- 
tions are averaged over the last 10 frames. 



We indicate by bo the initial mass present on solid bound- 
aries. The solid surface grows whenever b{r,t) exceeds the 
threshold value Z^max- The new solid node is chosen with uni- 
form probability among the surrounding nodes belonging to 
the liquid phase. Only direct neighbors are taken into account. 
The cumulated mass of the solid node that triggered surface 
growth is set to zero. 

Lattice Boltzmann simulations process and output data in 
model units. The different quantities will be expressed in 
terms of the basic model units for mass, length and time, 
denoted by mu, lu and ts, respectively. Direct comparison 
with physical systems is possible after suitable transforma- 
tions |22|. This procedure is not strictly necessary since, in 
order to describe properly the real systems, it is sufficient to 
have the same Reynolds number Re for fluid flow |39|, and 
the same Damkohler number Da and supersaturation xj/ for 
solute transport |fT9ll20ll22ll . These dimensionless parameters 
are the same in LB and physical units. For the reactive infil- 
tration of molten Si in a capillary of 10 [/im] in width, the 
Reynolds number can be estimated to be Re « 10~^ ifTolfTsl . 
This Reynolds number asks for extreme simulation condi- 
tions that can not be achieved. Indeed, a channel of width 
50 [lu] using = 1 [ts] would require a velocity of infiltra- 
tion M « 10^^ [lu/ts]. Thus, in the following, it should be 
borne in mind that our simulations significantly overestimate 
the typical Reynolds numbers for the problem at hand. In 
other words, this means that the process of infiltration is too 
fast. This numerical deficiency is partly reduced by the fact 
that the precipitation process from supersaturation should not 
have a strong dependence from the velocity of the invading 
front CSlEOlEa. 



III. RESULTS AND DISCUSSION 

A. Interfaces and surface tension 

Cohesive forces among molecules of the same species have 
the effect to give rise to interfaces separating the fluid com- 



ponents. Given two immiscible fluids, the work necessary in 
order to deform the interface of area A is 5W = ydA, where 
7 is the surface tension. By taking into account the work of 
pressure forces, the condition of mechanical equilibrium leads 
to Laplace law li40J 



AP=- in 2D. 
r 



(4) 



In the above formula, r is the radius of the spherical droplet. 
The pressure drop across the interface has the explicit form 
AP = Pin — Pout, with self-explanatory notation. In order to 
understand the interplay between interfacial tension, wetting 
and capillary flow, we first consider the formation of stable 
droplets for systems containing two fluid components. The 
size of the simulation domain is A^^ x A'^^ with A'^^ — 200 [lu]; 
periodic boundary conditions are applied. In the center is 
placed a droplet of fluid 1 immersed in fluid 2. Inside the 
bubble, pi is the main density and p2 is the dissolved density. 
For the initial condition, we assume that pi/pi = 2.5% and 
outside the bubble their values are inverted. The initial total 
density is denoted by po = pi +P2. By choosing the initial 
radius of the droplet to be r = N^/ it follows that both 
fluids have equal total masses, as assumed in Ref. BTI . The 
densities of both fluids determine the pressure P. Let pi and 
P2 be the densities of fluids 1 and 2 at point r — then 
the pressure is given by 125112611411 



P(r) 



Pi(r)+p2(r-) 



-Pi(r-)p2(r-) 



(5) 



In Eq. |4] for Pin we use the value obtained by averaging over 
a square of side length 2Lo + 1 centered in the simulation do- 
main; Lo = 20 [lu]. Instead, for Pout, the average is taken 
over the square [0,2Lo + 1] x [0,2Lo + 1]. Every dynamics 
amounts to lOO'OOO timesteps. For the analysis, 50 uniformly- 
spaced frames are collected. In order to derive the radius of 
the droplets, it is necessary to determine the location of the 
interface separating the fluid components. To do this, we start 
from the point {0,Ny/2) and consider at every step points in- 
cremented by one unit in x, i.e. x — ^ jc + 1 , applying every time 
the following rules. If pi(x+ l)/pi(x) > L05 holds, it fol- 
lows that the density for the first fluid increases at least of 5% 
in passing from x to x+l. We consider these lattice sites as 
belonging to the interface. For this set of points, we determine 
J^max and Xmin and it is assumed that the interface has coordi- 
nates X = (xmax +^min)/2 and y — Ny/2. Figure [T| illustrates 
the results according to the analysis proposed in Ref. [41 1. 
It turns out that the surface tension y is appreciable only for 
GcPo > 1, as thoroughly discussed in Ref. HTI . We consider 
outliers the two points for which G^Po < 1 and exhibiting 
a weak phase separation (broad interface), since inspection 
of the density values during the dynamics indicates that the 
systems are still evolving towards the equilibrium state after 
lOO'OOO timesteps. For Po = 2 [mu] and GcPo = 1 [lu/ts^], 
the average over the last 20' 000 timesteps for a simulation of 
500'000 timesteps yields pi = 1.17744 [mu]. For the shorter 
dynamics we obtained pi = 1.40671 [mu]. Furthermore, there 
appears that, even for the longer simulation, the density val- 
ues still exhibit systematic variations of at least 1% over the 
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Figure 2: Sessile droplet system for the determination of the contact 
angle 6. The solid line is the tangent at the triple line. Color code 
based on density. 



last 20' 000 timesteps. It is thus clear that Pi/po tends to- 
wards 0.5 141J. For I'OOO'OOO iterations this limit is further 
approached since pi = 1.12514 [mu]. On the other hand, for 
po 2 [mu] and GcPo = 1-4 [lu/ts^], we obtain pi = 1.83668 
[mu] in the case of a long simulation. The results remain prac- 
tically unchanged since for the shorter simulation we have 
Pi ~ 1.83733 [mu] and the difference between the two cor- 
responding values for the surface tension is on the order of 
10^^. For systems with higher surface tension the conver- 
gence towards the equilibrium state is faster. In the follow- 
ing, we will concentrate ourselves on three particular values 
of surface tension, leading in general to consistent and quite 
accurate results. 



B. Contact angle 

The wetting behavior of a surface against a liquid is usu- 
ally characterized via the equilibrium contact angle. Given the 
profile of the droplet, it is defined as the angle formed by the 
tangent at the contact line, that is, the edge where the phases 
solid (S), liquid (L) and vapor (V) meet (see Fig. |2]). For 
macroscopic, spherical droplets. Young equation [42] relates 
the contact angle 9 to the interfacial tensions Yij =S,L,V): 



cos 9 = 



Tsv - TsL 
7 



For the solid-liquid surface tension we continue to use the 
simpler notation 7. If < 90°, the liquid wets the substrate 
(hydrophilic regime); otherwise, if > 90°, the substrate is 
hydrophobic. 

In our simulations, the contact angle is determined using the 
sessile droplet method. The simulation domain is A'^v — 400 
[lu] long and Ny = 200 [lu] wide. The substrate coincides with 
the bottom of the simulation domain; the other boundaries are 
periodic. In the initial state, a square of side length 80 [lu] of 
the first fluid is centered in the simulation domain, placed in 
contact with the solid phase. The systems consist of a binary 
mixture where the droplet is represented by the wetting fluid 
immersed in the non- wetting component. As in Sec. Ill A 
for the initial condition, P1/P2 — 2.5% with po = Pi +P2 = 2 
[mu/lu^]. The simulations are performed with three different 
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Figure 3: Dependence of the contact angle 6 on the solid-liquid in- 
teraction parameter Gajs,! = ^Gajs ^ for ths simulation of the sessile 
droplet system (see Fig.|2j. Three different surface tensions are con- 
sidered by varying the parameter Gc (see Fig. [TJ. Points represent 
the results directly obtained from the profiles. The results predicted 
by Eq.|6] correspond to the dashed lines. The results of the simula- 
tions are averaged over the last 10 frames. The cases for which it is 
difficult to extract the contact angle are omitted. 



values for Gc. The equilibrium contact angle is varied by con- 
sidering different values for the interaction parameters Gads. a- 
The evolution of the systems consists of lOO'OOO timesteps. 

In the analysis, the contact angle is extracted using the 
method employed in Ref. B3l . In the micron scale, the droplet 
profile can be safely assumed to be spherical p4'|. Here, our 
systems refer to millimeter-sized droplets [9, 10 1. So, let h be 
the height of the spherical cap and R the base radius (contact 
line curvature), then the radius and the contact angle can be 
written respectively as 



2h 



and 9 



^ arccos 



(?) 



for the hydrophilic (R/h > 1) and hydrophobic (R/h < 1) 
regimes, respectively. If \ R/r \ > 1 we evaluate the function 
arccos using its Taylor series up to third order. For the de- 
termination of the height and base radius of the droplets the 
interface between the wetting and non-wetting fluid is derived 



as explained in Sec. Ill A An alternative way for calculating 



the contact angle is to use the formula BTll 



cos 9 = 



Gads,2 — Gads,] 



P1-P2 



(6) 



Pi and P2 are the main and dissolved densities averaged over 
the region [0,2Lo + 1] x [0,2Lo + Ijwith Lo = 5 [lu] centered 
in the point {Nx/2,h/2). Figure b] plots the contact angle 
versus the interaction parameter controlling adhesive forces 
between the solid and liquid phases. As the surface tension 
increases, it appears that the agreement between the two ap- 
proaches is more stringent. Our findings are in good accor- 
dance with those reported in Ref. Ii4ri . 



Figure 4: Set-up for the simulations of capillary-driven flow. Color code based on density: the first fluid is represented in red (high) and the 
second one in blue (low). Black points represent the solid boundary. The simulation domain has periodic boundary conditions. Top: Initial 
configuration. The first fluid fills the capillary along a distance of 50 [lu]. Middle: State of the system after the formation of the meniscus and 
partial penetration. The flat interface guarantees that the first fluid can be regarded as forming an infinitely large phase, mimicking in this way 
a reservoir. Bottom: In this case the color code is based on solute concentration. The region filled with the first fluid is at a higher concentration 
and the dividing interface from the low-concentration region approximately coincides with the interface for fluid flow. 
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Figure 5: Infiltration front at different times. Black points represent the profile as obtained from simulations. The green line is a fit to the data 
using a circle as model function. The red straight line is the tangent at the point where the contact angle 9 for the wetting fluid is determined. 
In general, the results determined in this way are in reasonable accordance with those based on capillary pressure (cf. Ref. L13J ). 
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Figure 6: Liquid displacement inside the capillary versus time. 
Points represent the simulation results while the solid lines corre- 
spond to the theoretical prediction of Eq.[8] The data shown in red are 
obtained with the surface tension y = 0.12338 [lu-mu/ts^] (Gc = 0.8 
[lu/mu/ts^]). For the results displayed in green the surface tension is 
7= 0.16403 [lu-mu/ts^] (Gc = 0.9 [lu/mu/ts^]). The results of the 
simulations are averaged over 10 frames. 



C. Capillary flow 

Capillary flow occurs when a liquid spontaneously spreads 
inside a channel. This phenomenon was exhaustively mod- 
eled by Washburn and Lukas 1451 l46l . It may be mentioned 
that historically capillarity played also an important role in 
the development of molecular theories 1471 . For two fluids 
with the same density and dynamic viscosity, by taking into 
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Figure 7: Dynamic contact angle, average value, as a function of the 
capillary width H (see Eq.[8](. The solid lines represent the equilib- 
rium contact angle as obtained using the sessile droplet system (see 
Sec ImB] ). 



account inertial effects, the centerline position of the invading 
front z{t) satisfies the differential equation Ii48^ 



d^z(0 12^ dz(f) 27COS0 



df2 //2p (Jf 



HpL 



(7) 



The first term on the left-hand side accounts for the inertial 
resistance and the second one for the viscous resistance. The 
term on the right-hand side is due to capillary forces. H and 
L are the height and length of the capillary. Q is the contact 
angle formed by the wetting fluid. Under the assumption of a 
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Figure 8: Typical interstice morphology as resulting from surface growth. The solid phase is represented in black. In red tonality the main 
component for mass transport and fluid flow. Top: Representation for solute concentration; Bottom: Representation for fluid density. 
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Figure 9: Time dependence of the infiltration depth with reactive boundaries. Data in red correspond to the surface tension 7= 0.12338 
[lu-mu/ts^] (Gc = 0.8 [lu/mu/ts^]). The data represented in green refer to the surface tension y = 0.16403 [lu-mu/ts^] (Gc = 0.9 [lu/mu/ts^]). 
The solid lines indicate the theoretical predictions in the absence of surface growth, Eq.[8] 



constant contact angle, a solution of Eq.|7]is obtained by ||T3]| 

Vcan/ZcOsS r / , \ , n 

z(t) = ?d [exp ( - f/fd) + f /fd -l]+zo. (8) 

The initial position of the interface is denoted by zo- Kap = 
7/ /i is the capillary speed and = pH^/12jj. is a typical tran- 
sient time 1 1 3 1 . 

For times sufficiently large, Eq. [sjleads to z(f) °^ t and the 
infiltration rate is given hy K = Vcap/Zcos 0/6L. Thus, this 
model reproduces correctly the linear dependence of the dis- 
placement front with time observed for molten Si and some 
alloys in porous carbon llO UTTI . Nevertheless, some hypoth- 
esis on which the result of Eq. [8] resides are not fulfilled in 
the experimental systems. For example, the liquid and vapor 
phases do not have approximately the same density. In the LB 
model, the density of both fluids determine the surface ten- 
sion (see Eqs.|4]and[5]). In the sequel, we will focus primarily 
on systems having an equilibrium contact angle around 30° 
(see Sec. 111B[ |. This result is consistent with the experimental 



results for molten Si and alloys on porous carbon substrates 
|P, TOl. Instead, the hypothesis assuming the same viscosity 
for both fluids is more restricting. Indeed, the attempt to in- 
crease the ratio between the liquid and vapor phases in the LB 
model would result in a breakdown of the linear Washburn law 
lfT2ll . This means that the present model only reproduces the 
macroscopic desired behavior and the effect of the reaction at 
the contact line is equivalent to the presence, in the capillary, 
of a second fluid as viscous as the wetting component. As a 
consequence, it follows that the role of the surface reaction 
between Si and C to form SiC can not be reduced to a sim- 
ple transition from a non-wetting to a wetting regime at the 
contact line of the invading front. 

In order to give evidence of the behavior predicted by 
Eq. [8] we consider systems similar to that shown in Fig. |4] 
As before, the initial densities satisfy pi/pi = 2.5% with 
Po = Pi + P2 = 2 [mu/lu^]. The width of the capillary varies 
from 25 to 100 [lu], corresponding to the height of the sim- 
ulation domain. The length of the simulation domain is in- 
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Figure 10: Az = Zt — Zr m the course of time. The lower script t refers to the theoretical value without reaction, Eq.[8] the symbol r indicates 
instead the front displacement for reactive boundaries. Red for data with surface tension 7 = 0.12338 [lu-mu/ts^] (Gc = 0.8 [lu/mu/ts^]); 
green for data with surface tension 7= 0.16403 [lu-mu/ts^] (Gc = 0.9 [lu/mu/ts^]). It should be noted that linearity indicates that pore closure 
occured in the case of reactive boundaries. This phenomenon is more frequent for narrow interstices. Another clear cause of linearity resides 
in the length of the capillary, that is, when the fluid reaches the other extremity in the absence of surface growth. 



stead set initially to I'OOO [lu]. The length of the capillary 
is L = 500 [lu]. The left extremity of the capillary coincides 
with X = 450 [lu]. Every system is let evolve for lOO'OOO 
timesteps. The dynamics is studied by collecting 50 evenly- 
spaced frames. In the analysis of the data, we first track the 
invading front along the centerline of the capillary, by em- 



ploying the method described in Sec. IIIB The contact angle 
6 in Eq.|8]is derived from the capillary pressure using the for- 
mula AP — 27COS 9/H 1121 [131 . In this case, the pressures Pin 
and fout are computed by averaging the pressure P{r), Eq.|5] 
over regions of the type [0,2Lo + 1] x [0,2Lo + !]■ These re- 
gions are centered around a point 20 [lu] apart from the cen- 
terline position of the meniscus in both the wetting and non- 
wetting fluids. It is assumed that Lq = 10 [lu]. This procedure 
yields the most accurate matching between theory and simu- 
lation results 1 13 1, allowing also to avoid difficulties related 
to the proper definition and determination of the profile near 
the contact line BTI . Figure |5] shows the invading front in the 
course of time for the system of Fig. |4] Agreement with the 
theoretical prediction is satisfactorily, as shown in Fig.|6] Of 
course, the surface tension remains unchanged since it is not 
a property depending on the geometry of the system. It turns 
out that the dynamic contact angle is on average larger than the 
equilibrium contact angle. Figure[7]indicates that the dynamic 
contact angle increases with H. Nevertheless, the process of 
infiltration is still faster for the capillaries of larger width (see 
Fig.|6]l. The present simulations are also repeated for L = 750 
[lu] (Af, = 1'500 [lu]) and L I'OOO [lu] (A^ = 2'000 [lu]) us- 



ing similar systems. The results are not presented for brevity. 
The analysis leads to the same conclusions discussed so far. 



D. Capillary infiltration witii reactive boundaries 

The SiC formation at the interface is responsible for wetting 
and spontaneous infiltration. This process is not accounted 
for directly since our systems are composed by a wetting and 
a non-wetting fluid, separated by an interface. This means 
that in the present modeling scheme the reaction at the triple 
line and related phenomena are not taken into account [10 17] 
i23J . Furthermore, solute transport does not affect fluid flow. 
Since the fluid and the solute represent the same substance, 
it follows that fluid flow is influenced only by surface growth 
and not by the diffusion and precipitation processes. It is thus 
assumed that the amount of Si involved in the surface reaction 
is negligible with respect to the part regarded as a liquid. 

The role of reactivity is investigated using again the set- 
up illustrated in Fig. |4] in the course of evolutions totaling 
lOO'OOO timesteps. Regarding fluid flow, we consider all the 
systems of Sec. IIIC| with the same initial conditions and 
the same settings defining the different wetting behaviors. 
For solute transport, in the initial condition, the region oc- 
cupied by the wetting fluid is filled with solute at the con- 
centration of Ci = 10^^ [mu/lu^]. In the region occupied by 
the non-wetting fluid, the solute concentration is set to C2 = 
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Figure 11: Time evolution of the maximal width Ae of the solid phase inside the capillary. Data in red have surface tension y = 0.12338 
[lu-mu/ts^] (Gc = 0.8 [lu/mu/ts^]); data in green have surface tension y= 0.16403 [lu-mu/ts^] (Gc = 0.9 [lu/mu/ts^]). 



2 • 10 ^ [mn/lu^]. For the parameters of Fs, Eq.js] we choose 
G, = -4.875 • 10-3 [mu/lu/ts2], = 1 and Co = 4.9 • 10"^ 
[mu/lu^]. The saturated concentration is quite small compared 
to the bulk value of solute: it is chosen to be Cs = 5 • 10^^ 
[mu/lu^], determining the relative saturation = Ci/Cs = 2. 
As a consequence, for relatively small values of ki-, the wetting 
fluid can be regarded as a reservoir for solute \2\. Further- 
more, these settings guarantee that in the non-wetting fluid, 
no solute deposits on the surface. The initial mass on solid 
boundaries is /jq = 2 • 10^^ [mu], and the threshold value for 
surface growth is assumed to be bj^ax = 10^^ [mu]. In order 
to adjust the effect of reaction relative to diffusion, it is suf- 
ficient to change the reaction-rate constant kj, without vary- 
ing the relaxation time for solute transport, set to Ts = 1 [ts] 
|fT9l |20l . Indeed, the systems can be classified according to 
the relative saturation = Ci/Cs and the Damkohler number 
Da = kyNy/D. The reaction-rate constant is varied accord- 
ing to the rule k, = {8/5)/T with T = 1'000,2'000, . . . ,5'000 
[ts]. Indicatively, this means that the flat surface start growing 
after 1'000,2'000, . . . ,5'000 timesteps. Figure [| shows a typ- 
ical configuration. It can be seen that the solid phase grows 
especially near the capillary throat. As evidenced by Fig. [6] 
the invading front varies linearly with time also for reactive 
boundaries. We see that by varying 7, k^ and L the relative 
position between the various curves appears to be the same. 
This means that the hierarchy between the associated infiltra- 
tion rates remains unchanged. It is of course necessary to take 



into account that narrow interstices obstruct sooner. In Fig. 10 



we compare more closely the position of the invading front 
to the theoretical prediction in the absence of surface reac- 



tion, Eq.[8] The retardation in terms of distance is obtained by 
considering Az — Zt — Zi-; the lower scripts t and r stand for the- 
oretical and reaction, specifying the two types of data. First 
of all, regarding the modeling approach, it is interesting to see 
that the curves corresponding to different surface tensions are 
separated more neatly. For the higher surface tension, more 
distance divides the invading front with reaction to that with- 
out reaction. Another general remark is that, for the wider and 
longer interstices, the retardation in terms of distance tends to 
increase faster and almost linearly towards the end of the sim- 
ulations. According to this representation of the simulation 
results, long interstices seem to be the optimal configuration, 
but this is not true. For example, for y— 0.12338 [lu-mu/ts^] 
and H —IQQ [lu], the invading front spans the distance of 400 
[lu] if L = 500 [lu] as opposed to approximately 300 [lu] if 
L — I'OOO [lu]. Let us now consider the maximal width Ae of 
the solid phase in the channel. This quantity can not exceed 
Ny/2. This value corresponds to pore closure. Figure 11 sub- 
stantiates the discussion in Ref. |fT9l related to the enect of 
the fluid velocity on boundary reactivity. It arises that these 
curves primarily depend on the reaction-rate constant k^. The 
surface tension 7 and the length of the capillary L do not seem 
to affect significantly the quantity Ae. Curves for interstices 
with different width H differ appreciably only when pore ob- 
struction occurs. As a result, the process of surface growth 
leading to pore closure appears to be independent of the in- 
filtration velocity. Narrow capillaries are thus more subjected 
to closure and their infiltration velocities are also smaller. It 
follows that narrow interstices are particularly detrimental for 
reactive fluid flow. 
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IV. CONCLUSIONS 

In the present work we have studied through LB simula- 
tions in 2D the capillary infiltration in a single interstice with 
reactive boundaries. The complex process of liquid Si infil- 
tration in C preforms was at the center of our attention. Se- 
vere inconsistencies exist between experimental and numer- 
ical conditions. For example, another simplifying assump- 
tion is that the LB simulations are performed in the isother- 
mal regime. It has been demonstrated that SiC formation is 
a highly exothermic reaction |8|. As a consequence, quanti- 
tative agreement is poor and it seems difficult to pinpoint the 
origin of inaccuracy given the mutual dependence of the vari- 
ous complex mechanisms at play. Nevertheless, the proposed 
analysis still provides us with some guidance. The thicken- 
ing behind the invading front effectively hinders the process 
of capillary infiltration. It turns out that wide and short in- 
terstices can limit the effect of pore closure (see Figs. [9] and 
[T0| . In other words, porous pathways as straight as possible 
should ease and dominate infiltration. Concerning the indus- 
trial manufacture of C/SiC composites (49], it follows that the 
most extended surface of the carbon matrix should be put in 
contact with molten Si, of course, if the porosity is isotropic. 



Then, surface growth does not result in a uniform corrugation 
of the initial flat surface but concentrates near the throat (see 
Fig. [8}. Thus, another way to reduce the slow-down of reac- 
tive capillary infiltration could consist in having round-shaped 
morphologies for the structure of the porous medium. Impor- 
tantly, narrow interstices are doubly disadvantageous since the 
process of pore closure can be regarded as independent of the 
infiltration velocity of fluid flow (see Fig. 11 1. In that respect. 



it is worth noticing that the benefits of bimodal distributions 
for pore size have been recognized |50|. Our next step would 
be to relate more thightly the coupled phenomena of the evo- 
lution of pore structure and capillary infiltration to the param- 
eters controlling the porosity. 
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